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ABSTRACT 

Cosmological constraints derived from galaxy clusters rely on accurate predictions 
of cluster observable properties, in which feedback from active galactic nuclei (AGN) is 
a critical component. In order to model the physical effects due to supermassive black 
holes (SMBH) on cosmological scales, subgrid modeling is required, and a variety of 
implementations have been developed in the literature. However, theoretical uncer- 
tainties due to model and parameter variations are not yet well understood, limiting 
the predictive power of simulations including AGN feedback. By performing a detailed 
parameter sensitivity study in a single cluster using several commonly-adopted AGN 
accretion and feedback models with FLASH, we quantify the model uncertainties in 
predictions of cluster integrated properties. We find that quantities that are more sen- 
sitive to gas density have larger uncertainties (~ 20% for M gas and a factor of ~ 2 for 
Lx at i?5oo), whereas Tx, Ysz, and Yx are more robust (~ 10 — 20% at Rsoo)- To 
make predictions beyond this level of accuracy would require more constraints on the 
most relevant parameters: the accretion model, mechanical heating efficiency, and size 
of feedback region. By studying the impact of AGN feedback on the scaling relations, 
we find that an anti-correlation exists between M gas and Tx, which is another reason 
why Ysz and Yx are excellent mass proxies. This anti-correlation also implies that 
AGN feedback is likely to be an important source of intrinsic scatter in the M gas -Tx 
and Lx~Tx relations. 

Key words: galaxies: clusters: general — hydrodynamics — methods: numerical — 
galaxies: active 



1 INTRODUCTION 

Clusters of galaxies are useful probes of cosmological param- 
eters, provided that their masses can be determined accu- 
rately from multi-wavelength observations calibrated using 
theoretical models. However, it is still a challenge for current 
theoretical models to reproduce all the observed properties 
of the baryonic content of clusters. Despite the fact that 
current cosmological simulations with radiative cooling and 
supernova feedback are able to reproduce profiles of the in- 
tracluster medium (ICM) outside the cores, the simulated 
cluster cores generally suffer from the over-cooling problem, 
e.g., the fraction of cool-core (CC) clusters and stellar frac- 
tion in these simulations are too high compared to observed 



values (see review by Borgani & Kravtsov ( 2009 ) and refer- 



ences therein). Therefore, some additional forms of heating 



Email: hsyang@umich.edu 



are required to suppress cooling. Feedback from active galac- 
tic nuclei (AGN) is one of the most promising candidates, 
as observations of X-ray cavities blown by jets from the cen- 
tral AGN imply a mechanical power that is comparable to 
the X-ray luminosity, suggesting a feedback loop might be 
at work ( |Dunn fc Fabian|2008| . 

Since a wide range of length scales is involved in this 
feedback loop, from the accretion disk of the supermas- 
sive black hole (SMBH) on parsec scales to clusters on 
Mpc scales, direct simulation with all relevant physics is be- 
yond current computational power. Thus cosmological sim- 
ulations with AGN feedback to date have had to model its 
sub-resolution effects by linking the resolvable scale (usu- 
ally ~kpc) to the SMBH accretion disk scale using some 
simplified assumptions ( Sijacki et al.|2007 Booth & Schaye 



2009||Dubois et al.|[20l6| |Gaspari et al.||2011[ >. These stud 
ies have had success in explaining the cosmic evolution of 
SMBH, star formation history, and local scaling relations 
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between the black hole mass and host properties, though 
current AGN models are still quite phenomenological. One 
important question to ask is whether these simulations can 
simultaneously reproduce observed ICM properties as well. 



2 METHODOLOGY 



2.1 Simulation Setup 



Analyses in this direction have started recently (Puchwein 
et al.|2008 Gaspari et al.pOlT l and need to be further ad- 
dressed. 

If cosmological simulations with AGN could accurately 
predict the core properties of clusters, they could provide 
crucial information for cluster cosmology, such as the CC 
fraction as a function of mass and redshift, or the parame- 
terizations of scaling relations. The CC fraction is important 
for understanding the X-ray selection bias toward CC clus- 
ters because of their peaked central surface brightness. For 
calibrating the scaling relations of multi-wavelength cluster 
surveys, one may not be comparing apples to apples if such a 
selection bias is not adequately accounted for. Moreover, the 
core properties may have an impact on the evolution in the 
slopes, normalizations, and scatter of the scaling relations. 
In order to obtain meaningful cosmological constraints using 



the self-calibration method (Levine et al. 2002 Majumdar 
fc Mohr|2004 Lima fc Hu|2004 l, the parametrizations of the 
scaling relations have to be informed by numerical simula- 
tions. 

However, the uncertainties in the existing AGN sub- 
grid models are not yet well understood. Since it is still 
unknown how to link the accretion and feedback across dif- 
ferent scales, there is a great amount of freedom to imple- 
ment and parametrize the AGN subgrid models. The model 
parameters sometimes do not have a clear connection to ob- 
servable quantities, and hence constraints from observations 
cannot be easily applied. Moreover, because these cosmolog- 
ical simulations require significant computational resources 
to run, it is difficult to perform detailed parameter studies to 
assess the robustness of the results. But in order to achieve 
predictions with high precision and controlled systematics, 
it is necessary to properly parametrize our ignorance in the 
AGN models and quantify the theoretical uncertainties. 

The aim of this study is to quantify the current theoret- 
ical uncertainties due to model variations in predicting the 
global ICM properties and thus provide a general guideline 
for cosmological simulations including AGN feedback. To 
this end, we implement a subgrid AGN unit in the FLASH 
simulation code that incorporates several existing AGN ac- 
cretion and feedback models. To study the effect of AGN 
feedback on cluster observables, we put these models in an 
idealized cluster and explore a wide range of parameters 
systematically. Connections between the model parameters 
and observable quantities are provided whenever possible. 
We identify the numerical details and parameters that the 
results are most sensitive to. For all the models that success- 
fully self-regulate black hole growth and reproduce observed 
cluster profiles, we then study the influence of AGN feedback 
on integrated cluster properties. 

The structure of this paper is as follows. The analytical 
and numerical approaches are described in § [2] The result of 
the sensitivity test is presented in § [3] In § [4J we first quan- 
tify the model uncertainties of integrated cluster properties, 
and then study the effects of AGN feedback on the scaling 
relations. Finally, we conclude and discuss possible ways to 
improve the subgrid models in § [5] 



et al. 20081. The cluster is set up in the same way as in 



We performed three-dimensional hydrodynamic simulations 
with radiative cooling and AGN feedback within an isolated 
cluster sitting in a 2048 kpc box using the adaptive mesh re- 
finement (A MR) code FLASH 3 ( |Fryxell et al.|2000| |Dubey 



Cattaneo & Teyssier (20071 and has properties similar to 



M87. The cluster gas is initialized assuming a polytropic 



equation of state (EOS) (Komatsu & Seljak 20011 and is 
in hydrostatic equilibrium in a fixed NFW ( jNavarro et al.| 



1995 1 gravitational potential. Self-gravity of the gas is not 
included. The cluster has a virial mass of 1.5 xlO 14 M©, con- 
centration of 5.53, and gas fraction of 0.1. This gas fraction 
is chosen to match the observed value for clusters of similar 
masses (e.g. Gonzalez et al. 2007). The contribution to the 
baryon fraction due to stellar mass is not included because 
it is small compared to the gas contribution. A 3 x 10 9 Mq 
black hole (BH) is placed in the center; it is only used for 
computing the accretion and feedback quantities and does 
not contribute to the gravitational potential. The base grid 
resolution is 32 kpc throughout the simulation volume, and 
the region surrounding the central black hole is refined pro- 
gressively to the maximum resolution (e.g. 1.0 kpc for the 
fiducial run). To accommodate the extent of feedback with 
high resolution, we define the maximally-refined region as a 
box of width 120 kpc for the bubble feedback model and 60 
kpc for the jet model (see 



2.2.2 for details of these mod- 



els). The diode boundary condition is used; this is similar to 
the outflow boundary condition but does not allow matter 
to flow into the domain. 

Radiative cooling is computed using Sutherland & Do- 
|pita| ( [T993| assuming 1/3 solar metallicity. Star formation 
and feedback from supernovae are neglected because they 
themselves have different implementations and require de- 
tailed comparisons of their own (e.g. |McCarthy et al.|20TTj ). 
In this study we intend to investigate the modeling of AGN 
alone and avoid confusion due to possible interference with 
other subgrid physics. We do not expect our main con- 
clusions to change because in most runs the gas densities 
never reach the conventional star formation threshold of 
jjh = 0.1 cm' 3 . Also, in our analysis of cluster integrated 
properties, we choose observables that are insensitive to the 
dense, cold gas surrounding the SMBH (see §[4|. Note also 
that in this paper we only focus on the hydrodynamic models 
employed in previous cosmological simulations, and hence 
the effects of magnetic fields are neglected. We refer read- 
ers to |Sutter et al.| ( |2012| ) for an investigation of magne- 
tized AGN feedback models. The Hubble constant h = 0.65 
is used. When overdensity quantities are quoted, they are 
computed using the overdensity radius Ra where the en- 
closed average density is A times the critical density of the 



2.2 AGN Subgrid Models 

Subgrid models in cosmological simulations have been de- 
veloped with varying levels of sophistication. Accretion rate 
calculations can vary from the simple Bo ndi (|Bondi||1952[) 
rate or its modifications ( Sijacki et al.|200T Booth & Schaye 
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20091, to accretion of cold gas only (Pizzolato & Sokc r|2005 
Gaspari et al.pOll l. To model quasars at high redshifts, it 
is commonly assumed that a fraction of their radiative en- 
ergy is transformed into thermal energy in the surround- 
ing gas (e.g. Di Matteo et al. 2005 Bhattacharya et al. 



2008). The mechanical input from AGN can be modeled 



using large-scale jets when the resolution permits (e.g. 


Cat- 


taneo & Teyssier 2007; Sternberg et al.||2007| |Dubois et al. 


2010 


Morsony et al.||2010| |Gaspari et al. 


2011 


). Since jets 



eventually inflate bubbles, it is easier computationally to 
place already-formed bubbles (e.g. Sijacki et al.|2007 Booth 
fc Schaye|[2009 l. 

Since it is infeasible for us to explore all the models pre- 
sented above, we select some representative accretion and 
feedback models. For estimating the BH accretion rate, we 
consider the a model proposed by [Sijacki et al.| ( |2007[ ). The 
feedback from the AGN is modeled based on the bubble feed- 



back of Sijacki et al. (20071 or the jet feedback of Cattaneo 
|fc Teyssierj ( |2007| ). Since our main goal is to quantify the 
model uncertainties, these models are chosen to cover very 
different methods of implementation and parameterization. 
We summarize the important aspects of these models in the 
following. 



2.2.1 SMBH Accretion 

To include the accretion onto the SMBH self-consistently 
in cosmological simulations, the simplest approach is to es- 
timate the accretion rate using the Bondi-Hoyle-Lyttleton 
I Bondi 1952) accretion rate: 



M BH oc Mb 



ndi = 47tG 2 m1hP/Cs, 



(1) 



where Mbh is the BH mass, and p and c s are the gas den- 
sity and sound speed, respectively. Cosmological simulations 
usually do not have sufficient resolution to resolve the Bondi 
radius, r Bon di = GM/c 2 a , as well as the multiphase gas when 
the density is high enough to trigger star formation. There- 
fore the density (temperature) at the Bondi radius would 
likely be higher (lower) than values on the grid. The actual 
Bondi accretion rate would thus be underestimated, which 
is reflected by the proportionality in the above equation. 

The a model assumes a constant coefficient, i.e., Mbh = 
«MBondi- Based on the resolution argument, ct = 1 is jus- 
tified for our simulated cluster because of its flat gas profile. 
However, the value of ct is often taken to be ~ 100 in previ- 



ous works to drive substantial black hole growth ( Di Matteo 
et al.|2005||Sijacki et al.|2007| [Bhattacharya et al.|2008| , in- 
dicating the large uncertainty in the estimation of accretion 
rates. In our parameter survey, we will vary a from 1 to 100 
in order to cover situations where the actual accretion rate 
is underestimated as well as where it is overestimated. 



1 It is possible to consider other forms of proportionality, such as 
the model proposed by Booth & Schayc (2009), which is consis- 
tent with the Bondi prediction when simulations have sufficient 
resolution or when gas densities are lower than the star formation 
threshold nj/ =0.1 cm -3 ; otherwise the proportionality is den- 
sity dependent to account for accretion of multiphase gas. Since 
in our simulations the gas densities seldom reach the star forma- 
tion threshold, the /3 model is equivalent to a = 1 in our current 
setup. 



After computing the accretion rate, an upper limit is 
imposed corresponding to the Eddington rate, MEdd = 
(47rGA/BHTrip)/(efOTc), where m p is the mass of the proton, 
or is the Thompson cross-section and et is the radiative 
efficiency. 

The region for computing the accretion rate and the 
region from which to deplete the accreted gas are typically 
set to a few zones around the central BH, but their sizes 
are essentially arbitrary. We denote these two radii as J? aC c 
and Rdep, and we probe several values in § |3.2| Note that 
during strong accretion events, gas in a grid cell can possibly 
be completely removed and cause an unphysical surge of 
gas around the black hole. In such cases, we increase the 
depletion radius by a small amount so that no more than 
10% of the gas on a grid cell is removed in one timestep. 
We expect this condition to have negligible effects because 
it does not occur frequently, and the results are insensitive 
to the depletion radius (see 



3.2 l 



2.2.2 AGN Feedback 

The feedback from the AGN to the surrounding gas is then 
computed according to the accretion rate. There has been 
growing observational evidence for an anti-correlation be- 



tween radio loudness and SMBH accretion rate (Ho 2002 



Sikora et al. 20071. That is, radio jets are associated with 
systems having lower accretion rates, while objects with 
higher accretion rates, like quasars at higher redshifts, are 
radiatively efficient, analogous to states of X-ray binaries 
( Fender et aL]|1999| |Gallo et aL]|2003[). For t his reason, we 
follow the prescription in Sijacki et al. (20071 for switching 



to the quasar mode when the accretion reaches 1% of the 
Eddington rate. In the quasar mode, the radiative energy is 
thermally coupled to the surrounding gas: the energy depo- 
sition rate is E = e r e{MBHC 2 , where e r is the quasar heating 
efficiency. The region into which to dump the quasar thermal 
energy is chosen to be four zones in radius, though it is arbi- 



trary. Note that in some other subgrid models (e.g. Booth & 
Schayc 2009 ) there is no division between the quasar and me- 



chanical feedback since they are both assumed to be purely 
thermal and spherically distributed. However, we will show 
that when feedback energy is injected in the thermal form, 
the size of the region has a significant impact on the results. 

At low accretion rates, one can choose either bubble or 
jet feedback. For the bubble feedback ( |Sijacki et al.||200"7| ), 
the feedback energy is distributed in terms of thermal energy 
within a spherical region around the SMBH. Bubbles are 
only formed when the BH mass increases by a fraction 5bh 
since the last bubble formation. When a bubble is formed, 
purely thermal energy is injected: 



E = e m £fC 8 Mbh, 



(2) 



where e m is the efficiency of mechanical heating, and 
5Mbh = <5bhMbh is the increase in BH mass since the last 
bubble was formed. The injected energy is distributed in a 
mass-weighted sense within a sphere of radius 

1/5 

R = R 



E/E 
p/Po 



(3) 



where the scaling parameter values _Ro, Eo, and po are mo- 
tivated by observed bubble sizes. The bubble centers are 
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randomly displaced within a sphere of radius i?dis centered 
on the black hole. In the fiducial run Rdi s = R- We also ex- 
periment with cases where bubbles are fixed at the central 
black hole, i.e., Rdis = 0. 

In contrast to bubble feedback, which only injects ther- 
mal energy into the surrounding ICM, the jet feedback sim- 
ulations inject some or all energy in kinetic form. The jet 
models are not intended to simulate the relativistic jet di- 
rectly, but rather the non-relativistic outflow from the ac- 
cretion disk |Proga 2007) or decelerated large-scale jet af- 
ter it has entrained some intergalactic medium during its 
propagation ( |Feretti et al.||1999| [Eaing fc Bridle|2002[ |. The 



Table 1. Survey of Numerical Parameters in the Bubble Model. 



ratio of injected thermal to kinetic energy depends on the 
parametrization and is different from study to study. While 



Gaspari et al. (20111 



and Dubois et al. ( 


20101 


|Cattaneo & Teyssier 


(|2007| 



jected energy is mostly thermal, depending on the amount 



of mass loading. The model of Cattaneo & Teyssier (20071 
is motivated by the observation that more massive, slow 
jets should couple more thermal energy with the surround- 
ings as they propagate. However, a purely kinetic jet in their 
model would produce relativistic velocities, which cannot be 
treated adequately in non-relativistic hydrodynamic simula- 
tions. Therefore, we modified their model and used a more 
general parametrization to allow tuning of the thermal to 
kinetic ratio, while the jet velocities are independently de- 
termined by the amount of mass loading. 

Our generalized jet model can be summarized as follows. 
The injection rates of mass, momentum, and energy onto 
the grid are treated as source terms in the hydrodynamic 
equations and are calculated by 



M = 7?Mbh|*|, 
P = v^yyer (1 - e m )A/ BH ctf, 



(4) 



E 



EfMaHC 1*1 



where r\ is the mass loading factor and E is the sum of 
injected thermal and kinetic energy, i.e., E = E t ^ + E^. For 
non-relativistic jets, Ey_ = P 2 /2M; thus in this model e m 
of the total energy goes into the thermal energy, and the 
remainder is kinetic. The jet velocity is c-y/2ef (1 — e m )/?7, 
which is ~ 10 4 km/s for i) = 100, et = 0.1 and e m = 0. The 
function ^ determines the spatial extent of the jet: 



*(x) 



2 < 



exp 



2 i 2 

x +y 



(5) 



The jet is aligned with the 2-axis, a nd the fe edback is applied 



to 



regions with \z\ ^ h cq and \J x 2 + y 2 



^ r c j. We also 



normalize the window function \& so that the total injected 
energy in the cylinder sums up to E. Note that there is 
no threshold for injecting the jets, so the jet feedback is 
continuous rather than episodic. 

Note that our jet model includes several modifications 
to that in |Cattaneo fc Teyssier| ( [2007[ ). In addition to changes 
in the parametrization as described above, we normalized 
the function as in their subsequent papers based on the 
same model ( |Dubois et al.|20To| ) . They fixed the mass of the 
black hole for computing the accretion rate, while our black 
hole can grow from gas accretion. This has negligible effects 
since the accretion rate is small in their setup. Also, we allow 
gas depletion and use a smaller radius for computing the 



Name Ax (kpc) 



TOTT 
BIB 
B1C 
BID 
B1E 

Scaling 
~E2K 



Varying Resolution 

03 r 

1.0 1 

2.0 1 

4.0 1 

8.0 1 



Alpha with Resolution 
2T0 2 



Table 2. Survey of Numerical Parameters in the Jet Model. 



Name 


Ax (kpc) 


h c j r c j (kpc) Rs.cc R dcp 


(zones 


Varying Resolution 


J1A 


0.5 


2.0 2.5 2 


2 


JIB 


1.0 


2.0 2.5 2 


2 


J1C 


2.0 


2.0 2.5 2 


2 




Scaling Jet Size with Resolution 




J2A 


1.0 


4.0 5.0 2 


2 


J2B 


2.0 


8.0 10.0 2 


2 


J2C 


4.0 


16.0 20.0 2 


2 


J2D 


1.0 


8.0 10.0 2 


2 




Varying 


Radii for Accretion and Depiction 




J3A 


1.0 


2.0 2.5 1 


1 


J3B 


1.0 


2.0 2.5 4 


4 


J3C 


1.0 


2.0 2.5 2 






accretion rate. When accounting for the above differences, 
we are able to reproduce their results. 



2.3 Model and Parameter Variations 

Since the details of each subgrid model are so different that 
it is infeasible to explore every implementation and param- 
eter, in this study we focus on those aspects of these models 
which are least constrained. To this end, the bubble and 
jet feedback models are chosen because they are different in 
many aspects, including the form of injected energy (thermal 
vs. kinetic), shape of injection region (spherical vs. jet-like), 
and periodicity of feedback (episodic vs. continuous). Com- 
paring these two distinct models allows us to understand the 
extent of current theoretical uncertainties due to these AGN 
models. 

We first explore 'numerical' parameters that are re- 
quired in the numerical implementations but are essentially 
arbitrary and often chosen based on numerical rather than 
physical considerations. Variations of these parameters are 
listed in Table [l] and Table [2] for the bubble and jet model, 
respectively. The tests are divided into groups, each of which 
examines one of the numerical parameters. We first do a con- 
vergence test by varying the peak resolution, Ax, for both 
the bubble model (group Bl) and the jet model (group Jl). 
Because the constant a model is motivated by the argument 
that the accretion rate is under-estimated due to resolution 
insufficient to resolve the Bondi radius, in group B2 we scale 
the a parameter with resolution to see if the results are con- 
sistent. In the jet model we do not vary a because we will 
show that the accretion rate is suppressed from the begin- 
ning and hence varying a has a minor effect. Since in the jet 
model the feedback is distributed in the very inner few kpc 
around the black hole, it is more likely to interfere with the 
peak resolution and the choice of the size of the region for 
gas accretion and depletion. Therefore, we experiment with 
scaling the jet size with resolution in group J2 and varying 
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Table 3. Survey of Physical Parameters. 



Name 


a 


Feedback 




&H(%) 




Region 


Varying Accretion 


P1A 


1 


bubble 


0.1 0.2 


0.01 


Ro 


= 30 h- L kpc, R dis = R 


P1B 


10 


bubble 


0.1 0.2 


0.01 


Ro 


= 30 h~ 1 kpc, R dis = R 


PIC 


100 


bubble 


0.1 0.2 


0.01 


Ro 


= 30 h- 1 kpc, R dis = R 








Varying Mechanical Heating 


Efficiency 






P2A 


1 


bubble 


0.1 0.02 


0.01 


Ro 


= 30 h~ L kpc, R diB = R 


P2B 


1 


bubble 


0.1 0.5 


0.01 


Ro 


= 30 h- 1 kpc, R dis = R 








Varying Feedback Frequency 






P3A 


1 


bubble 


0.1 0.2 


0.001 


Ro 


= 30 h~ L kpc, R diB = R 


P3B 


1 


bubble 


0.1 0.2 


0.1 


Ro 


= 30 h^kpc, R dis = R 








Varying Size of Feedback Region 






P4A 


1 


bubble 


0.1 0.2 


0.01 


Ro = 


= 15 h ^kpc, Rdi s = R 


P4B 


1 


bubble 


0.1 0.2 


0.01 


Ro 


= 5 h _1 kpc, Rdis = R 


P4C 


1 


bubble 


0.1 0.2 


0.01 


Ro 


= 15 h ikpc, Rdi s = 








Varying Thermal to Kinetic Ratio 






PSA 


1 


jet 


0.02 







= 2.5 kpc, h c j = 2.0 kpc 


P5B 


1 


jet 


0.02 0.5 





r «i 


= 2.5 kpc, h c j = 2.0 kpc 


P5C 


1 


jet 


0.02 1 







= 2.5 kpc, h c j = 2.0 kpc 


P5D 


1 


jet 


0.1 







= 2.5 kpc, h ei = 2.0 kpc 


P5E 


1 


bubble 


0.1 0.2 


0.0001 


R = 


= 2 h 1 kpc, Rdis = 



the accretion and depletion radii in group J3. The impacts 
of these numerical considerations will be discussed in § |3.2| 
Note that for these jet runs the feedback and thermal effi- 
ciency parameters are chosen to be ef =0.1 and e m = 0% 
(purely kinetic), as commonly adopted in previous jet mod- 
els ( |Gaspari et al.|2011| |Dubois et al.||2010 l. 

Table |3| summarizes the model or parameter variations 
that are more physically motivated. Again they are sepa- 
rated into groups. The effects we investigate here include 
varying the accretion strength, the mechanical feedback ef- 
ficiency, the threshold for triggering a feedback event, the 
size and center of the feedback region, and the ratio be- 
tween thermal and kinetic energy. Here we note that despite 
the difference in the parameterizations used in the bubble 
and jet models, we do expect their results to overlap when 
bubbles are injected almost continuously (with a small <5bh) 
into a small region centered on the black hole (run P5E), 
and when the jet is purely thermal (run P5C). The differ- 
ences between these two runs would mainly be due to the 
shape of the injection region. We discuss the influence of 
these 'physical' parameters in § |3.3| 

When not specifically mentioned, the parameters are 
fixed to their fiducial values: ef = 0.1, e r = 0.05, e m = 0.2, 
Ro = 30 h,- 1 kpc, E = 10 55 erg, p = 10 4 h 2 M© kpc" 3 , 
r c j = 3.2 kpc, and h e j = 2.5 kpc. Note that these values 
are chosen to match those adopted in previous cosmological 
simulations ( |Sijacki et al.|[2007 Dubois et al.pOlO i, which 
are able to successfully reproduce the observed black hole 
density, star formation history, and Lx~Tx relation. There- 
fore it is also one of our aims to examine whether the same 
set of parameters could also recover the properties of the 
ICM within observational constraints. 



3 SENSITIVITY STUDY 

In this section we present a sensitivity study of all the rel- 
evant parameters in the AGN subgrid models. Since most 
of these parameters are not well constrained due to lack 
of knowledge of the detailed physical processes, we will first 
start from a relatively large parameter space and then exam- 



ine whether the results are consistent with observed limits. 
By performing the sensitivity study we would like to iden- 
tify those variables to which the ICM properties are most 
sensitive. At the same time, this study also provides us with 
information about which ICM properties are most robust to 
the uncertainties in the AGN subgrid models. 



3.1 The Fiducial Run 

Here we show the general features of run PI A as an instruc- 
tive example. This run uses the bubble model with param- 
eters that are commonly adopted in previous cosmological 
simulations ( Sijacki et al.|2007 1. The left column of Figure [l] 
shows the evolution of the SMBH accretion rate, the mass of 
the black hole, the power of feedback (bubble energy divided 
by the duty cycle), and the duty cycle (time between two 
feedback events). The right column shows the profiles of gas 



density, temperature, pressure and entropy (K = T/rhe ) at 
different times. For this run, the accretion rate starts from a 
value of ~ 6 x 10" 4 M yr _1 (M B a/M Edd ~ 9 x 10~ 6 ). The 
black hole grows slowly in its mass and generates weak bub- 
bles only every few hundred Myr for the first ~ 4 Gyr. The 
injected energy delays the strong cooling which would occur 
if there were no feedback. But it is still unable to balance ra- 
diative losses, so the accretion rate grows rapidly and reaches 
1% of the Eddington rate, triggering quasar-mode feedback 
at t ~ 5 Gyr. At this time the accretion rate is held near the 
value 1 Mq yr _1 by powerful and frequent (~ 0.1 — 1 Myr) 
feedback events (both bubble and quasar modes), which in 
turn heat and expand the surrounding gas and cause the ac- 
cretion rate to drop. The corresponding decrease of feedback 
then results in another round of strong cooling and feedback 
events. The cluster then fluctuates with a timescale of ~ 3 
Gyr till the end of the simulation. 

The gas profiles change according to the AGN activity. 
During the first ~ 4 Gyr, the central density increases and 
the temperature decreases due to radiative cooling, similar 
to the case without AGN feedback. But once powerful feed- 
back from the AGN starts to take place after t ~ 5 Gyr, the 
gas entropy within the injected bubbles is raised, so that 
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AGN outbursts, e.g., the large floor and entropy inversions 
at t = 6 Gyr and t = 12 Gyr. These powerful events are also 
reflected in the density profiles, which have more flattened 
cores than observed. Note that the bubbles generated in this 
run have radii ~ 100 — 200 kpc, larger than typical sizes of 
observed X-ray cavities (McNamara & Nulscn 20071. There- 
fore, the influence of the bubbles may be overestimated due 
to the parameter choice of the bubble size (i.e. Ro in Eq. 
[3f. As we will discuss in § |3.3.4| though this problem can 



be alleviated by choosing a smaller Ro, it means that the 
default parameters for bubble sizes may not be applicable 
to every cluster but need to be fine-tuned. 

This run shows that such AGN feedback models can 
successfully self-regulate black hole growth, as also demon- 



strated in previous work (e.g. Sijacki et al. 20071. It also 
produces cluster profiles that are in general consistent with 
observations. In the following subsections we will start vary- 
ing the parameters to see how they affect the evolution of the 
AGN activities and cluster profiles. Note that our simula- 
tions do not include gas self-gravity. Though fragmentation 
of gas due to self-gravity could be important during the hi- 
erarchical formation of clusters ( Dubois et al.|2012 l, for the 
already-formed cluster in our simulations, the effect of gas 
clumping is minor because the Jeans length is large, except 
in the central regions where catastrophic cooling occurs for a 



few runs (see § 3.3.5 1. Also, including self-gravity would pos- 
sibly deepen and flatten the cluster potential well and result 
in a flatter asymptotic pressure profile. However, since the 
main interest of this study is the differences due to model 
variations, and gas only contributes to a minor portion of 
the total cluster potential, we do not expect the inclusion of 
gas self-gravity to significantly change the results. 



Figure 1. Results for the fiducial run P1A. Left column (from 
top to bottom): Evolution of black hole accretion rate, black hole 
mass, power of feedback, and duty cycle. The X-ray luminosity 
inside the core (R ^ 0.15-Rsoo) is overplotted with the injected 
power using the dashed line. Right: Radial profiles of gas density, 
temperature, pressure, and entropy. Grey areas are the observed 
ranges of density, pressure and entropy profiles fro m|Croston et al.| 
| |2008[ l, [Arnaud et al.| | |2010[ | and |Cavagnolo et al.| l |2009| ), respec- 
tively. 



the gas is heated and pushed outward. After later times the 
cluster oscillates around a quasi-static profile that is flatter 
and hotter than the initial profile. 

For the density, pressure and entropy profiles, we over- 
plot our results with the observed profiles recently compiled 



for a large sample of clusters by Croston et al. ( 2008 1 , Ar- 



|naud et ah] ( |2010j ) and |Cavagnolo et al.| ( |2009| ), respectively. 
These observed profiles are remarkably uniform and self- 
similar at outer radii, while the dispersion increases toward 
the center. We note that despite these energetic AGN out- 
bursts, the pressure profiles at all times lie well within the 
observed range. The universality of the pressure profiles is 
maintained because the density and temperature of the bub- 
bles compensate each other to reach hydrostatic equilibrium 
with the surrounding ICM. 

On the other hand, the entropy profiles, though follow- 
ing a standard 'power-law plus floor' profile in general, some- 
times contrast with the observed profiles right after powerful 



3.2 Numerical Parameters 

In this section we present results for varying the numerical 
parameters both in the bubble model (Table [T]) and the jet 
model (Table For the bubble model, we study the ef- 
fect of peak resolution, as well as scaling the a parameter 
with resolution. For the jet model, we test variations of the 
resolution, jet sizes, and the radius for computing accretion 
rate and for removing gas. The definitions and explanations 
for these parameters can be found in § |2.3| We probe the 
sensitivity to these parameters by examining the evolution 
of the black hole accretion rate. The cluster profiles are not 
plotted here since they are closely related to the accretion 
rate, as seen in the previous section. 

The convergence test for the bubble model is shown 
in Figure [5] (top left). We allow a large range of variation 
for the peak resolutions (as large as 8 kpc) because current 
high-resolution cosmological simulations can typically reach 
resolutions of a few kpc, but to go beyond that is progres- 
sively more difficult due to computational costs. Therefore it 
is important to understand whether such simulations with 
subgrid AGN models are numerically converged. We find 
that as the peak resolution is increased, fluctuations with 
shorter timescales and with larger amplitudes are captured, 
whereas for runs with degraded resolutions the accretion 
rates react more slowly, reach the first peak later, and have 
less variation. This makes the SMBH in the lower-resolution 
runs increase its mass at a later time, but it grows to a larger 
value at the end of the simulation by time 12 Gyr. In par- 
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Figure 2. Black hole accretion rates for varying numerical parameters (see Table [T] and Table 121. Top left: Varying the peak resolution 
(group Bl) and scaled a with the resolution (group B2) for the bubble feedback model. Top right: Affecting the accretion rate by varying 
the peak resolution in the jet feedback model (group Jl). Bottom left: Effects of different jet sizes with fixed/scaled peak resolution 
(group J2). Bottom right: Results of changing the radius for computing the accretion rate and the radius for removing the accreted gas 
in units of zones (group J3). 



ticular, the final black hole mass for run B1E (Aa; = 8 kpc) 
differs from run BIB (Aa; = 1 kpc) by about a factor of 2. 
Therefore simulations with poorer resolutions may under- 
estimate the variations in accretion rate and corresponding 
ICM properties. They may also find more massive black hole 
populations when other conditions are held the same. 

Since the constant a in the bubble model is often in- 
voked to compensate for the underestimation of the accre- 
tion rate due to resolution, in principle a larger a should 
be used when the resolution is poorer. In group B2 we test 
whether such scaling of a would help account for the dif- 
ference in resolution. As shown in Figure [5] (top left), the 
initial accretion rate for run B2A (Aa; = 2 kpc, a = 2) is 
boosted compared to the run with the same peak resolution 
(B1C). However, since this cluster has a flat initial entropy 
profile (recall that M B ondi oc n c /cf oc n c /T 3/2 oc K~ 3/2 ), 
this boost is not necessary to match the accretion rates of 
the higher-resolution runs (B1A and BIB). This points out 
one problem with the a accretion model, which is that a 
single constant value of a may not be appropriate for a pop- 
ulation of clusters with various core profiles. At later times 
when cooling and feedback events get stronger, the evolu- 
tion becomes nonlinear and depends sensitively upon the 
detailed interactions between the bubbles and quasars with 
the surroundings. So the outcome of run B2A is close to nei- 



ther run B1C nor run BIB. Therefore, simply scaling a by a 
constant factor does not in general work to compensate for 
the change in resolution. In fact, we will see in §[3X1] that 
choosing the value of a is nontrivial and has a great impact 
on the evolution of SMBH and ICM properties. 

For the jet feedback model, we also find more variation 
in the accretion rate for runs with higher peak resolution 
(Figure|5J top right). Despite the difference in the amplitude 
of fluctuations, the mean accretion rates are more robust to 
the resolution than in the bubble model, at least when the 
zone size is larger than 1 kpc (run JIB and J1C). 

The dramatic change in the accretion rate for run J1A 
may be understood in combination with the results of scaling 
jet sizes with resolution (Figure [2] bottom left). This prob- 
lem of large amplitude fluctuations occurs in all the runs 
where the jet size is much greater than the size of a grid 
cell (J2A-J2D). Recall that the radius for computing the 
accretion rate is set to 2 zones. Therefore when the accre- 
tion radius is small compared to the region used to apply jet 
feedback, such accretion rate estimates are sensitive more to 
the details in the feedback itself than the actual accretion 
from the ICM. The only exception is run J2C, which has 
very extended jets that have essentially the same effects as 
the bubble feedback. In particular, this run produces results 
that are comparable to bubbles with high frequencies (run 
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P3A) and small sizes (run P4C) that will be shown in later 
sections. 

Finally we vary the radius for computing the accretion 
rate, J? a cc, and the radius for removing the accreted gas, 
Rdcp (Figure [5] bottom right). We find that either changing 
the accretion and depletion radii, or removing gas depletion 
altogether, has a minor effect on the results. Note however 
that run J3A crashed at t ~ 8 Gyr with a sudden drop in the 
accretion rate. Since in this run the accretion radius is only 
one zone and is smaller than the jet size, the accretion rate is 
very sensitive to the central few zones within the feedback 
region. A slight displacement of gas within the innermost 
cells causes a sudden reduction in the accretion rate and the 
associated feedback, which induces an unphysical surge of 
gas into the central zones. Therefore, for numerical stability 
we recommend using an accretion radius larger than the 
region of jet feedback. 

To summarize briefly, we find that increasing the peak 
resolution generally produces more variable accretion rates. 
Bubble feedback suffers a greater influence when varying the 
resolution, whereas jet feedback is more robust, as long as 
the accretion radius is carefully chosen (larger than the jet 
size). 



3.3 Physical Parameters 

3.3.1 Dependence on Accretion Models 

Table [3] lists the variations of the physically-motivated pa- 
rameters under consideration. For the first group of runs we 
vary the method of computing the accretion rate. As can 
be seen from Figure |3j the evolution of the cluster is very 
different for different values of the accretion strength pa- 
rameter, a. For run PfA where a = 1, the accretion starts 
from a small value and thus the initial feedback power is 
small compared to the core X-ray luminosity. Here strong 
cooling occurs and triggers cycles of feedback events at later 
times, as described in § |3.1| On the other hand, in run PIC 
(a = fOO), the feedback power in the beginning is already 
comparable to the X-ray luminosity, so the cluster never goes 
through the strong cooling phase and is roughly in equilib- 
rium throughout the simulation time. The cluster profiles 
respond to the AGN activity in the same way as in the fidu- 
cial run, which is the reason why the gas properties have 
more fluctuations in run P1A than Pf C. 

The evolution of the cluster core can be further quanti- 
fied using the entropy floor ( Cavagnolo et al.|2009 \ , slopes of 
the density profile ( Croston et al.|2008 l and entropy profile 
(Sanderson et al. 20091, or the cooling time (Mittal et al 



2008}. Following the definition in |Mittal et al.| 02008} , we 
categorize cluster cores as strong cool cores (SCC; t coo i < 1 
Gyr), weak cool cores (WCC; f Gyr < t coot < 7.7 Gyr), 
and non cool cores (NCC; t CO oi > 7.7 Gyr). Based on this 
definition, the cluster starts with a WCC with t coo \ ~ 7.5 
Gyr. In all the bubble models studied here, the cluster never 
reaches the SCC state because the overall heating is very ef- 
fective. For run P1A, the cooling time drops to its minimum 
of tcool ~ 2 Gyr at t ~ 4 Gyr and climbs up to a NCC in 
the end. Run PIB is similar to P1A. For run PIC, the cool- 
ing time instead keeps rising, so the cluster became a NCC 
cluster soon after the start of the simulation. 

The large influence of the assumed accretion model on 



the evolution history of the cluster poses a great concern for 
simulations including AGN feedback. It implies that given 
the uncertainties in the accretion mechanisms, current AGN 
subgrid models have very limited power to predict the evo- 
lution of cluster core properties. That is, simulations with 
different accretion models can produce very different results, 
e.g., the fraction of CC versus NCC clusters as a function of 
time. Note that the cosmological simulation of |Sijacki et al.| 
(20071 uses a = 100, which is very effective in suppressing 



the formation of cool cores as we have shown. Indeed their 
simulations generally overpredict the fraction of NCC clus- 
ters at the present day compared to observations]^] There- 
fore, in order to produce robust results, accurate modeling 
of the accretion onto the SMBH is crucial. 



3.3.2 Varying the Mechanical Heating Efficiency 

The mechanical heating efficiency, e m , parametrizes how 
much feedback energy is actually converted into thermal en- 
ergy and used as a source of heating. Observationally the ra- 
tio of cavity power to the Bondi accretion rate is estimated 



to be a few percent (Allen et al. 20061, which motivates 



previous simulations to adopt similar values for the net effi- 
ciency, £{e m . However, it is still unclear how the cavity power 
is converted into heat and how long this process takes. If the 
bubbles do not mix with the ICM efficiently, in principle 
the mechanical heating efficiency could be much lower. For 
example, Vernaleo & Reynolds ( |2007| l used purely hydro- 
dynamic simulations and estimated the fraction of injected 
kinetic energy going into internal energy (i.e. e m ) to be only 
a few percent. But since the actual magnitude of mechanical 
heating would depend on the details of mixing, simulations 
with more realistic physical treatment of the ICM are re- 
quired to pin down this number. Here we probe the range 
e m = 0.02 — 0.5, which is permitted by current constraints 
and covers values commonly used in previous simulations. 

Figure [4] shows the SMBH and cluster evolution with 
varying e m (group P2). In general the changes in e m do not 
alter the fate of the cluster, in contrast to the variation in the 
accretion models discussed in the previous section. For all 
three runs, the cluster goes through gradual cooling for the 
first 4 Gyr, which eventually grows and triggers a sequence 
of feedback events, just like the fiducial run P1A. The more 
powerful bubbles for run P2A (e m = 0.5) only delay the time 
of strong cooling a little bit, but for the first 4 Gyr the results 
are almost indistinguishable. After 5 Gyr, the cluster again 
starts to oscillate among states that are closely related to the 
feedback activity. When the mechanical heating efficiency 
is larger, the feedback is less frequent, more powerful and 
more effective in reducing the accretion rate. The reduced 
accretion rate takes a longer time to grow back to 1% of the 
Eddington rate, so the cluster fluctuates with a timescale 
of ~ 5 Gyr for run P2A (e m = 0.5), while for run P2B 
(e m = 0.02) the cluster oscillates with smaller amplitude 
and timescale. 

Despite the overall similarity of behavior, the change 
in e m does result in several noticeable trends. First, the 
growth of the black hole mass is sensitive to the efficiency. 
For larger efficiencies, the accretion rate is suppressed so 
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that the black hole does not grow as much as when the ef- 
ficiencies are smaller. Second, efficient mechanical heating 
produces large density plateaus and entropy floors in the 
cores. Particularly at moments right after an energetic out- 
burst, the density and entropy profiles can have temporary 
excursions that go beyond observed ranges. This problem 
may be alleviated by reducing the scaling of the bubble ra- 



dius (see § 3.3.4|. Lastly, for runs with small efficiencies, the 



temperature profiles sometimes have a peak at the center. 
This is due to the concentrated thermal energy injection in 
the quasar mode, which is invoked more frequently in these 
cases to help the bubbles at high accretion rates. This is an 
issue shared with all the models which input thermal energy 
into a small region, such as those jet models with nonzero 



thermal components (see § 3.3.51. To avoid this undesirable 
feature of too much quasar-mode feedback, we did a test run 
identical to P2B (e m = 0.02), but with the switch to quasar 
mode turned off (i.e., the bubble mode is applied through- 
out the simulation). We find that pure bubble heating alone 
with such small efficiencies is insufficient to halt the cooling 
catastrophe. The accretion rate increases to ~ fOO Mq yr -1 
at t ~ 5 Gyr and generates a huge bubble that essentially 
blows all the cluster gas away. This puts a lower limit on 
the mechanical feedback efficiency, as also found by |Gaspari| 
et all (12011 1. 



3.3.3 Varying the Feedback Frequency 



In the bubble model employed in this paper (Sijacki et al 



|2007j), the bubbles are inflated when the black hole mass 
increases its mass by a fraction <5bh- Therefore, a larger 5bh 
corresponds to a longer duration between successive bubble 
events, as can be seen from Figure [5] The exact relation- 
ship between the parameter 5bh and the duty cycle actually 
can be derived given the criterion for triggering bubbles. For 
this particular bubble model, the increase in BH mass be- 
tween subsequent AGN outbursts is JbhA/bh — MbhAt. 
Therefore, 



Mbh 
10 9 Mr. 



M BH 



1 M yr- 1 



Myr. 



(6) 



Indeed, we can see from Figure [5] that the initial duty 
cycle scales with Sbh , and that the evolution of the duty cy- 
cle is roughly inversely proportional to the accretion rate. 
For M87, if we take observationally constra ined values, 



0.026 M yr" 1 , Mb h ~ 3 x 10" M 
^ yr (e.g. |Million et al.|2010 



10' 



(Allen et al. 



we obtain 



M BH 

|2006[ ), and At 

Sbh ~ 0.01%, which is consistent with the value used in the 
fiducial run. Note that Eq. [6] is only true for this specific 
bubble model. The duty cycle would have different scalings 
with BH mass and accretion rate when a different crite- 
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Figure 4. Evolution of AGN activity and cluster profiles for different mechanical heating efficiency t n 
P1A). 



(group P2 plus the fiducial run 



rion is used to trigger bubbles. For example, in the model 
of Booth & Schaye (20091, the injected energy per bubble 



is fixed according to a minimum heating temperature, so 
that E — constant oc Mbhc 2 At gives At oc Mgg with no 
explicit dependence on BH mass. On the other hand, the 
model of |Battaglia et aL (20101 sets Ar = constant instead. 



Therefore, how to trigger bubbles in the AGN models is not 
completely arbitrary, but in principle can be constrained by 



observed scalings of duty cycles (e.g. Shabala et al. 20081, 
though these measurements are themselves very difficult. 



Fortunately, we find that the frequency of injections in 
the model is not critical for the results. As shown in Fig- 
ure [5] the amplitudes of variation in the accretion rate and 
gas properties are similar for different <5bh- The increase 
in temperature and entropy is slightly higher for smaller 
Sbh, mainly because the injected energy is distributed in a 
smaller region, as the bubble sizes scale with the injected 
energy, which is smaller for smaller <5bh (Eq.(2j. But com- 
pared with the effects of accretion and mechanical heating 
efficiency, varying Sbh, or the frequency of bubble injection, 
does not have as large an impact on the overall evolution of 
SMBH and the ICM. 



3.3.4 Dependence on the Region of Feedback 

Here we explore the effect of varying the region of bubble in- 
jection (group P4), including the scaling coefficient for bub- 
ble radii, Ro, as well as the displacement from the central 
AGN, i?dis- As in the fiducial case, run P4A and P4B in- 
ject bubbles whose centers are randomly displaced within 
a sphere of radius J?dis, but with smaller bubble radii. The 
typical size of bubbles is 100-200 kpc, 50-100 kpc and 20-30 
kpc for the fiducial run PI A, run P4A and run P4B, respec- 
tively. As shown in Figure [5] reducing the bubble sizes is 
very effective in suppressing the accretion rate because of 
concentrated heating. The influence on the evolution of the 
cluster core is even more than increasing the heating effi- 
ciency (Figure |4|. In other words, the evolution of SMBH 
and ICM properties is very sensitive to the choice of bubble 
radii. 

Compared with the fiducial run, the cluster profiles for 
run P4A (_Ro = 15 h~ l kpc) have a smaller flattened core 
since the energy injection is more concentrated. Note that 
now the density and entropy profiles at all times are consis- 
tent with the observed range. This implies that the choice 
of bubble radii is not completely arbitrary. Bubbles that are 
too large may produce density and entropy profiles that are 
inconsistent with observations. 
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Figure 5. Dependence on the threshold for bubble injection (group P3 plus the fiducial run P1A). The threshold, <5bh> ' s defined to be 
the minimal fractional increase in black hole mass required to inflate a bubble. Smaller thresholds generate more frequent bubbles. 



The sizes of the bubbles cannot be too small either. 
For run P4B (Rq — 5 hT 1 kpc) where small bubbles are 
randomly injected around the black hole, though concen- 
trated heating greatly slows down cooling and accretion, 
catastrophic cooling still occurs at t ~ 10 Gyr and gen- 
erates bubbles that dramatically heat the core. This may be 
due to the fact that the radius within which the cooling time 
is less than 10 Gyr is around 100 kpc. So the bubbles only 
heat a small fraction of gas within the cooling radius, while 
a substantial amount of gas is still allowed to cool and flow 
to the center at later times. 

We also did one run with the center of bubbles fixed on 
the central AGN (run P4C, Rdis ~ 0). Other parameters are 
the same as in run P4A (7?o = 15 /i _1 kpc, Ra ls = R). From 
Figure [6] we can see that these two runs produce almost 
identical results, except that fixed bubbles tend to produce 
less smooth profiles than randomly-positioned bubbles. But 
these differences are minor. Therefore, as long as the size of 
bubbles and thus the input energy density are comparable, 
where around the SMBH to dump the energy has a lesser 
effect. 

In summary, the size of the region used to inject thermal 
energy (but not so much the displacement from the black 
hole) is crucial in predicting the evolution of the SMBH and 



the ICM. Bubbles that are too large push too much gas 
outward and raise the entropy floor to an unrealistic level, 
while bubbles that are too small may not be able to heat all 
the region where it is needed, allowing catastrophic cooling. 
Therefore, for the current bubble model, or any model that 
requires setting the size of energy injection by hand, it would 
be difficult to find one parameter or scaling that works for all 
clusters. And even if the results are permitted by observed 
limits, the predictions would still be very sensitive to the 
chosen bubble sizes. 



3.3.5 Effect of Thermal to Kinetic Ratio 

In this section we explore the models where the feedback 
energy is discharged in the form of jet, which has a shape 
function (Eq. [5]l aligned with the z-axis, as opposed to the 
spherical bubbles discussed in previous sections. In our gen- 
eralized parametrization of the jets (Eq. |4]| , the amount of 
thermal energy and kinetic energy can be tuned using two 
parameters: ef, the feedback efficiency, or the ratio between 
total injected energy and the rest mass energy of the SMBH, 
and e m , the fraction that goes into thermal energy. The com- 
parison of different thermal to kinetic ratios (group P5A- 
P5D) is displayed in Figure [7] 
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Figure 6. Effect of the region for bubble injection (group P4). Rq is the scaling coefficient for bubble radii as in Eq. [3] i?dis is the 
displacement of bubble center from the AGN. The evolution of the SMBH for the fiducial case P1A is plotted using black curves, but its 
cluster profiles are omitted here (see previous figures). 



The first thing to note for the jet models is that the ac- 
cretion rates are immediately reduced to ~ 10~ 4 M© yr -1 
as soon as the simulation starts and are more so when the 
thermal efficiency e m is greater. This suppression is due to 
the fact that in the jet models, energy is injected within 
only a few kpc around the AGN, instead of large bubbles 
that extend up to tens or hundreds of kpc. As shown in the 
previous section, decreasing the size of the region used to dis- 
tribute thermal energy can suppress the accretion rate very 
effectively. Therefore, the jets just resemble tiny bubbles. In 
order to verify whether the bubble and jet models are con- 
sistent under similar conditions, we did a run where bubbles 
with fixed radii 2 ft -1 kpc are generated almost continuously 
(run P5E). Figure [8] shows that this run indeed reproduces 
the case P5C of purely thermal jets of identical net feed- 
back efficiencies, efe m = 0.02. The small fluctuations for run 
P5E just reflect the fact that bubbles are produced every 
few timesteps rather than perfectly continuously. This test 
demonstrates that the bubble and jet models are numerically 
consistent and are degenerate when appropriate parameters 
are chosen. Moreover, it again emphasizes the point that the 
choice of the size for energy injection in the AGN subgrid 
models is nontrivial. 



Runs P5A-P5C compare different ratios of thermal to 
kinetic energy, with the total feedback efficiency ef kept fixed 
at 0.02. The fraction of energy that goes into thermal energy 
is 0%, 50%, and 100% for run P5A, P5B, and P5C, respec- 
tively. As shown in Figure [7] all these jet models in general 
produce similar overall evolution of the accretion rates and 
cluster profiles. The accretion is halted by concentrated feed- 
back from the beginning so the black hole grows very slowly 
for the first few Gyr. Since the injected power is smaller 
than the X-ray luminosity, the cluster gas cools more and 
more rapidly to a cool-core state at t ~ 6 Gyr, The rapid 
cooling quickly feeds the black hole and increases its mass, 
which allows the jets to stabilize cooling afterwards. The 
only exception is run P5A, which fails to overcome catas- 
trophic cooling at t ~ 10 Gyr. So self-regulation of black 
hole growth may not be achieved by purely kinetic feedback 
with efficiencies that are too small. 

As expected, the accretion rate is initially more sup- 
pressed for higher thermal efficiencies. But interestingly, af- 
ter t ~ 6 Gyr run P5B (e m — 50%) becomes the most effec- 
tive. In other words, the most effective way to stifle cooling 
is not necessarily to dump all the feedback energy in ther- 
mal form, but rather a combination of thermal and kinetic 
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Figure 7. Varying the thermal and kinetic efficiencies in the jet feedback model (group P5A-P5D). tf is the feedback efficiency, i.e., the 
ratio between total injected energy and the rest mass energy of the SMBH, and e m is the fraction that goes into thermal energy. 



feedback that facilitates mixing of the heated gas with the 
surroundings. 

We also performed another run P5D with purely kinetic 
feedback but higher total feedback efficiency, £f = 0.1. Like 
run P5A, which is also kinetic but with tf = 0.02, the initial 
accretion rate is not affected in the beginning since the ki- 
netic energy has not transformed into heat. Their differences 
become more evident as feedback energy is thermalized and 
as the jets become more powerful after t ~ 6 Gyr. Note that 
at later times, the level of suppression is comparable to run 
P5B. So both raising the total feedback efficiency and tun- 
ing the thermal to kinetic ratio can slow down the accretion 
and black hole growth. 

It is also instructive to compare run P5D (ef = 0.1, 
purely kinetic) with run P5C (ef = 0.02, purely thermal). If 
their results are comparable, it would imply that 20% of the 
kinetic energy is converted into heat, or a mechanical heating 
efficiency of 0.2, which is around values commonly assumed 
in AGN subgrid models using purely thermal feedback (| Si-| 
|jacki et aL 2007 Booth & Schayc 2009]) . Recent cosmological 
simulations also have found that either using purely kinetic 
feedback or purely thermal feedback assuming 15% for the 
mechanical heating could match the local black hole prop- 
erties ( Dubois et al.pbll l. Here we find that run P5D has 



somewhat better ability than run P5C to halt cooling, im- 
plying possibly a mechanical heating efficiency higher than 
20%. Note that small discrepancies are expected because of 
different sizes of thermal feedback used by different simula- 
tions. For our jets the feedback region is confined within the 
small shape function, whereas simulations mentioned above 



used either extended bubbles (Sijacki et al. 12007 1 or the 
neare st SPH particle ( |Booth fe Schaye||200"9j|Dubois et al.| 
20111. Note also that these arguments are based on simple 



hydrodynamic simulations, and would change if gas mixing 
is modified by additional physical effects, such as subgrid 
turbulence, viscosity, etc. 

The resulting cluster profiles again reflect the ability of 
jets to stop cooling. Since run P5B is the most effective, its 
temperature decreases and density increases the least. How- 
ever, in general the jets do not expel the gas or heat the 
gas as much as the bubbles, maintaining the core in the CC 
state. Note that the density profiles show a dense core of size 
~ 10 kpc at later times. This is due to a torus of cold gas 
that forms around the black hole when cooling is happening 
rapidly. In reality this cold gas should keep condensing to 
unresolved scales and form stars. This unphysical accumu- 
lation of cold gas is treated in previous work using different 
methods, including removing the cold gas with a sink term 
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4 6 
Time(Gyr) 

Figure 8. Comparing the bubble and jet feedback models. Run 
P1A is the fiducial bubble case, which generates large, randomly- 
displaced bubbles. Run P5C and P5E have the same net feedback 
efficiencies, £fe m = 0.02. But run P5C uses purely thermal jets, 
while run P5E uses tiny, fixed, and almost continuous bubbles to 
mimic the jet run P5C. This figure demonstrates that the bubble 
and jet models are numerically consistent when appropriate pa- 
rameters are chosen. Moreover, their differences are mainly caused 
by the size of the feedback region. 



in the continuity equation ( jGaspari et al.||2011[ ), or using 
an effective equation of state appropriate for multiphase gas 



( Dubois et al.||2010p . Since our simulation does not include 
these treatments, we will avoid deriving quantities that are 
sensitive to the central densities. For example, we will only 
use core-excised instead of total X-ray luminosity in later 
sections when we discuss cluster observables. 

Another thing to note is that any jet with nonzero ther- 
mal efficiency (run P5B and P5C) would produce a hot spot 
surrounding the black hole and hence a peak in the temper- 
ature profile near the center, which is also found by|Gaspari| 



et al. (20111 and Dubois et al. (20111. Therefore, we ad- 



vise simulators using concentrated thermal feedback such as 
thermal jets or quasar feedback to be cautious about this 
numerical effect when interpreting results near the region of 
feedback. 



4 IMPLICATIONS FOR CLUSTER 
OBSERVABLES 

4.1 Robustness of Integrated Properties 

As seen in the previous section, specific models or param- 
eters chosen can result in quite discrepant predictions for 
the evolution of the cluster profiles. Although the influ- 
ence of AGN feedback is strongest in the core region, we 
may ask whether, under the influence of this feedback, the 
global cluster properties can still preserve the observed sear- 
ings. The robustness of integrated quantities (e.g. measured 
within i?5oo) is particularly crucial for cluster cosmology. 
Since current constraints are often derived using calibra- 
tions of the mass-observable relations informed by numer- 
ical simulations, it is necessary to quantify the systematic 
uncertainties due to incomplete knowledge of the details of 
AGN feedback processes. 

The first question we wish to address is whether any of 



the model variations explored in the previous section predict 
global cluster properties that violate the observed scaling 
relations. To this end we compute several observable quan- 
tities integrated within a sphere with radius -R500, includ- 
ing the gas mass M s , spectral-like temperature T s i (Maz- 
zotta et al.||2004[ ), X -ray luminosity Lx, integrated Comp- 



ton y parameter due to the Sunyaev-Zel'dovich (SZ; Sun 



yaev & Zeldovich ( 1972 l) effect Ysz, and its X-ray analog 
Y x = M S T X ( |Kravtsov et al.||2006[ ). Since the spectral-like 
temperature and the X-ray luminosity are very sensitive to 
dense and cold gas, we excised the core region (< 0.15i?5oo) 
in order to avoid numerical effects due to the cold gas accu- 
mulated around the SMBH as described in § |3.3.5| For run 
P4C {Ro = 5 h,- 1 kpc) and run P5A (e f = 0.02, e m = 0%), 
the evolution after 9 Gyr is not included because these runs 
encounter the cooling catastrophe. 

Figure [9] compares the trajectories of observables to the 
scaling relations for the model variations explored in § |3.3| 
Each column compares results for a particular group in Ta- 
ble [3] including variations in the accretion model, mechan- 
ical heating efficiency, feedback frequency, region of feed- 
back, and thermal to kinetic ratio. From top to bottom we 
show the Mg-T s i, Lx~T B \, Ysz-Lx, and Ysz~Yx relations 
and overplot the observed relations and scatter for the REX- 
CESS samp le ( |Croston et al.|2008| [Pratt et al.|2009||Arnaud 
et al.|[2010 |. We note that the offsets in the normalizations 



are due to the different ways to derive the intrinsic and 
observed quantities. Firstly, the observed values of -R500 in 
these studies are obtained by matching the empirical M500- 
Yx relation in |Arnaud et al. (20071. As these authors point 
out, their hydrodynamic mass JV/500 may underestimate the 
true mass. Therefore it is likely that their -R500 is smaller 
than that used in our computation, which lowers the values 
of M g and Lx, but increases T s i, in a direction that could 
explain the shift in normalizations. Moreover, the X-ray lu- 
minosity in our calculation is bolometric and hence would be 
higher than the observed values, which are integrated over 
the energy range 0.1 — 2.4 keV. More detailed simulated ob- 
servations are required for direct comparisons between the 
simulations and observations. 

We find that despite the variation in the predicted clus- 
ter profiles produced by different subgrid models, the inte- 
grated properties for all the models evolve with amplitudes 
that are consistent with the observed scatter. In other words, 
when cooling is regulated, all the subgrid models are able to 
preserve global cluster properties as observed. This result 
gives us some confidence in the AGN subgrid models em- 
ployed in cosmological simulations. However, it also implies 
that these various models and parameters cannot be dis- 
tinguished by constraints on the integrated quantities of a 
single cluster alone, but must be constrained by observa- 
tions that are more sensitive to cluster cores, or by compar- 
ing scaling relations of a sample of simulated clusters with 
observations. 

Although none of the individual trajectories in Figure 
[9] violates the observed scaling relations, when compared 
against each other, the predictions for a particular observ- 
able at a particular time can still vary significantly among 
models. Thus the second question to ask is: how large are the 
theoretical uncertainties due to different AGN subgrid mod- 
els and variations in their parameters? Since AGN feedback 
is expected to have more impact at smaller radii, we compute 
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Figure 9. Trajectories of integrated observable properties on the scaling relations for variations of physical parameters explored in § |3.3| 
The rows from top to bottom show the M s -T s \, Lx-T s \, Ysz~Lx, and Ysz—Yx relations, respectively (see text for detailed definitions). 
The columns from left to right plot the runs with varied accretion model (group PI, Figurepl, mechanical heating efficiency (group P2, 
Figure 4 I, frequency of feedback (group P3, Figure[5]l, region of feedback (group P4, Figure |6|, and thermal to kinetic ratio (group P5, 
Figure 7 1, respectively. Overplotted are observed relations (solid) and r.m.s. scatter (dashed) for the REXCESS sample (i Croston et al.| 
|2008||Pratt et aI.|2009||Arnaud et al.|2010| . This fi gure illustrates that despite the discrepancies in the detailed evolution among different 
subgrid models as seen in previous figures, the integrated cluster observables still evolve with amplitudes that are consistent with the 
observed scatter. See Figure [12] and [13] for the detailed evolution of observables. 



the observables measured within several commonly-quoted 
overdensities, including -R2500, ^1000, ^500, and -R200, and 
compare their values with the fiducial runs, i.e., run P1A 
for the bubble model and run P4D for the jet model. The 
relative dispersion for a given observable O is computed by 
AO = \0— Ofiduciai I /Oflduciai ■ The results for five observables 
(T s i, M g , Lx, Ysz, and Yx) with varied subgrid models are 
presented in Figure [l"0| 

Comparing the five observables, the X-ray luminosity 
has the largest uncertainties due to variations in subgrid 
models (note that its plotting range is — 200%), and the 
gas mass is second. The other three variables (T s i, Ysz, and 
Yx) are more robust. When different groups of model varia- 
tions are contrasted, we find that the mechanical heating ef- 
ficiency and the size of the feedback region cause the largest 



variations in the predicted observables. The influence of the 
accretion models and feedback frequency is smaller, and the 
thermal to kinetic ratio has the least impact. 

As expected, since feedback from the AGN is more in- 
fluential toward the central SMBH, the model uncertainties 
are biggest for observables measured within i?2500- When 
quantities are integrated out to -R200, the uncertainties be- 
come small for all observables except the X-ray luminosity, 
because a large fraction of the total luminosity still comes 
from regions near the core. In Figure [TT1 and Table [| we 
summarize the maximum uncertainties among all models 
for each cluster observable versus the overdensity radius. 
We find that in general observables that are sensitive to gas 
densities are more poorly predicted. The total gas mass pre- 
dicted by different models has uncertainties ranging from a 
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Figure 10. Uncertainties of cluster integrated properties due to AGN subgrid model variations (see Figure [9] for explanations of each 
column). Plotted are the predictions of T s \, M g , Lx, Y$z, and Ky (top to bottom; the core is excised for computing T s ; and Lx ; see 
text for details) relative to the fiducial runs as a function of four overdensity radii within which the observables are integrated, including 
^?2500i KiouOi R500 and -f?200- Each color corresponds to a specific run, and each line represents the result at a given simulation time. 
Note that the plotting range for the X-ray luminosity is — 200%, since it has the largest uncertainty due to model variations. T s i, Ygz, 
and Yx are more robust. Comparing different groups of model variations (i.e. by columns), we find that the mechanical heating efficiency 
and size of feedback are the most influential, while feedback frequency and thermal to kinetic ratio play a minor role. 



few percent at R200, to ~ 20% at -R500, to ~ 100% at i?2soo- 
Thus the X-ray luminosity, which is proportional to density 
squared, can vary by factors of a few for all radii. The level 
of uncertainties is smaller and comparable for T B \, Ysz, and 
Y x , ranging from ~ 40 - 50% at i? 25 oo, to ~ 10 - 20% at 
#500, to ~ 5 — 10% at 7?2oo- 



4.2 Impact of AGN feedback on the Scaling 
Relations 

Keeping in mind the model uncertainties of integrated prop- 
erties shown in the previous section due to different evolu- 
tion in each model, next we study some general trends pre- 



dicted by all the models. In particular, we probe the impact 
of AGN outbursts on the cluster observables using cross- 
correlations among them. Note that in this study we focus 
on the global observable properties, that is, whether there 
will be observable features beyond the core due to the dis- 
turbances introduced by the central AGN. The influence on 
the core properties by AGN has been discussed extensively 



in previous work (see McNamara & Nulsen (20071 and ref- 



erences therein) and will be a part of our future work. 

Figure 12 (top row) plots the evolution of the AGN 
power, X-ray luminosity, and spectral-like temperature for 
all the bubble runs listed in Table [3] The jet runs are not 
shown here because their AGN power cannot be compared 
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Figure 11. Model uncertainties for cluster observables as a func- 
tion of overdensity radii. Values correspond to the maximum val- 
ues across all models shown in Figure [To] 



Table 4. Model Uncertainties for Cluster Observables Measured 
within Various Overdensities. 
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directly in the thermal form. However, their results can be 
well represented by the jet-like bubble run P5E, as discussed 
in section § |3.3.5| (see also Figure [8]). For all the models, the 
black hole self-regulates its growth when its feedback power 
is sufficient to balance the radiative losses by the cluster. 
However, the feedback power fluctuates around the mean 
after t ~ 6 Gyr with different amplitudes depending on the 
initial configurations and growth at earlier times. These fluc- 
tuations are present in the AGN activity as well as in the 
cluster observables. Taking run P2B (red curve) as an illus- 
tration, the strong AGN outbursts at t ~ 5 — 6 Gyr raise 
the entropy of the gas and hence induce the decrease in 
luminosity and increase in temperature at t ~ 6 — 7 Gyr. 
Similar effects can be seen for the peak in AGN activity 
around t ~ 10 Gyr and corresponding fluctuations in the 
observables at t ~ 11 Gyr. 

We therefore compare the AGN power with the changes 
in luminosity and temperature after t = 6 Gyr (bottom row 
in Figure [l2| to see whether their fluctuations are correlated. 
The notations dlog(Lx) and dlog(T s i) represent logarithmic 
deviations from their initial values; the feedback power is 
plotted with respect to the initial luminosity too. The cor- 
relation coefficients are given by the Spearman Rank-Order 
Correlation test (Press et al. 111992 1. As expected, there is a 



negative (positive) correlation between the feedback power 
and the time derivative of the luminosity (temperature). We 
note that the correlations are much weaker if the AGN power 
is compared with the instantaneous luminosity and temper- 
ature, because of the phase shifts between the peaks of AGN 



outbursts and the delayed responses of the ICM. Neverthe- 
less, since the luminosity and temperature react to the feed- 
back in phase, they have a strong anti-correlation as the 
system moves on the Lx~T a \ plane, as shown in the lower 
right panel in Figure |l2"| 

Recall that the ranges of trajectories on the Lx~T B \ 
plane predicted by all the AGN subgrid models are compa- 
rable to the observed scatter (Figure [9j. This implies that 
AGN feedback can drive a significant amount of the observed 
scatter in the Lx—T B \ relation, because of the anti-correlation 
between luminosity and temperature during the feedback 
events. This is in contrast to other physical processes such 
as cluster mergers, which tend to move the clusters along 
the scaling relations ( Yang et al.|2010 i. 

Similarly, an anti-correlation exists between M g and T s \ 



(Figure 13 left panel), which may also contribute to the 
scatter in the M s -T a \ relation. Moreover, the anti-correlation 
implies that Ysz and Yx , which are essentially the products 
of M g and T s i, do not deviate from the mean scaling relations 
significantly during AGN outbursts. As shown in the right 
panel of Figure (13} the variation in both Y parameters from 
the mid point is dlog(V) ~ 7.5%. Given the observed slope 



of the Y-M relation of 1.82 ( |Arnaud et al.||2007| l, it cor- 
responds to an implied uncertainty in mass dlogM ~ 4.1% 
if the Y-M relation is used to get a cluster mass. This is 
smaller than the mass uncertainties inferred from other ob- 
servables estimated by the same means, indicating that the 
Y parameters are robust mass proxies even under the strong 
influence of energetic AGN outbursts. 

Since we have demonstrated that the scatter in the Lx~ 
T B \ relation can be induced by feedback events, in Figure 
right panel) we correlate the Lx~T B i scatter with the 
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AGN power. Since we do not have a sample of clusters to 
derive the mean scaling relation, the scatter is computed by 
taking the logarithmic deviation from the observed relation 
shown in Figu re [9| As expected from the correlations found 
earlier (Figure 12 1, a negative correlation exists between the 
scatter and the AGN power. However, again the trend is not 
prominent because of the phase shifts. 

This result may have implications for observational 
studies that attempt to connect the Lx-T a \ scatter to the 



AGN radio power. Croston et al. (20051 found that radio 



loud AGN preferentially lie below the Lx~Tx relation as 
evidence for AGN heating. However, a more recent study by 



|Jetha et al. ( 2007 1 found a weaker relation. For illustration 
we plot the epochs when the AGN power is 0.5 dex more 
(less) than the zero-point value in filled (open) symbols (for 
clarity only data points at multiples of one Gyr are shown). 
As can be seen in the left panel of Figure [l4"l there is no clear 
segregation on the Lx~T s \ plane between the more powerful 
and the more quiescent populations, because the correlation 
between the scatter and AGN power is not strong enough 
(right panel). If the radio loudness is (roughly) proportional 
to the power of AGN, this may explain why observationally 
it is difficult to find a strong correspondence for a sample of 
clusters. 

The scatter in the Lx~Tx relation has long been known 
to be dominated by the core properties of clusters; exclud- 
ing the emission inside the core region can significantly re- 
duce the scatter (e.g. Pratt et al.|2009 |. Moreover, CC clus- 
ters generally have a higher normalization on the plane than 
NCC clusters. Here we explicitly show in Figure [l5]that such 
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Figure 12. Top: Evolution of the AGN power, X-ray luminosity, and spectral-like temperature measured within [0.15-l]i?5oo- All the 
bubble runs in Table [J are shown in different colors {black: P1A, gray: P1B, gold: PIC, blue: P2A, red: P2B, navy: P3A, cyan: P3B, 
green: P4A, pink: P4B, magenta: P4C, yellow: P5E). Bottom: Correlations for data points in the time interval 6 ^ t ^ 12 Gyr (r is the 
correlation coefficient) between the AGN power and the change in X-ray luminosity (left), the AGN power and the change in spectral-like 
temperature (middle), and the X-ray luminosity and spectral-like temperature (right). These correlations are expected because AGN 
outbursts result in reduction of the luminosity and heating of the cluster with slight time delays. 
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would depend upon model selection. Interestingly, cosmolog- 
ical simulations including bubble feedback have shown that 
AGN feedback is capable of suppressing the Lx~T s i nor- 
malization for low-mass clusters and steepening the slope 
to match observations ( jPuchwein et al.|2008[ ). However, our 
study suggests that there can be systematically different re- 
sults if one chooses different models, such as models with 
smaller mechanical heating efficiency, or the jet models. 



Figure 13. Trajectories on the M g — T s \ (left) and the Ygz~Yx 
(right) relations for 6 ^ t ^ 12 Gyr for the same runs as in 
Figure [l2| (r is the correlation coefficient). Similar to the Lx~ T s \ 
relation, M g and T s \ are anti-correlated. As a result, there is a 
tight positive correlation between Ygz an d Yx- 



a trend can be caused by the effects of AGN. During the 
AGN feedback events, the Lx~T B \ scatter is anti-correlated 
with the cooling time (right panel). Thus the CC and WCC 
clusters tend to lie on the upper half of the relation com- 
pared to NCC clusters (left panel). Since the luminosity and 
temperature studied here are core-excluded, even stronger 
trends are expected to be found for the core-included Lx~ 
T B \ relation. 

Note however that the NCC clusters at later times in 
the simulations are mostly produced by models in which 
the AGN feedback is either very powerful (high mechani- 
cal heating efficiency) or very extended (large bubble sizes), 
whereas the CC clusters are produced by the jet models. 
Therefore the exact amplitude of this segregation of CC and 
NCC clusters, or the suppression of Lx~T s i normalization, 



5 DISCUSSION AND CONCLUSIONS 

Feedback from the AGN is a crucial ingredient in modeling 
the observable properties of galaxy clusters. In the literature 
there has been a variety of AGN subgrid models employed 
in cosmological simulations. However, systematic parameter 
surveys and comparisons among different implementations 
are critical for understanding the robustness of their pre- 
dictions. In this study, we implemented several commonly- 
adopted accretion and feedback models into FLASH and 
systematically explored various parameters in an idealized 
cluster atmosphere. We first performed a sensitivity test of 
these subgrid models using a spectrum of parameters in or- 
der to understand their relative importance. We then quanti- 
fied the theoretical uncertainties of cluster integrated prop- 
erties due to model variations, and studied the impact of 
AGN feedback on the scaling relations by summarizing the 
results among different models. 

Since it is infeasible to explore every implementation 
and parameter of existing AGN subgrid models, our study 
is focused on two common approaches: to inject thermal en- 
ergy in an extended region to mimic already inflated bubbles 
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Figure 14. Trajectories on the Lx~ T & \ relation (left) and correlation between its log scatter and the power of AGN in the time interval 
6 t ^ 12 Gyr for the same runs as in Figure [l2] (right; r is the correlation coefficient). Outbursts that are 0.5 dex more (less) powerful 
than the mean are marked in red filled circles (blue open circles). More powerful AGN preferentially have smaller scatter (lie below the 
mean); however, the correlation is diluted by the phase shift between the AGN power and the observables shown in Figure [12| 
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Figure 15. Trajectories on the Lx~T s \ relation (left) and correlation of its long scatter with the cooling time in the cluster core in 
the interval 6 ^ t ^ 12 Gyr for runs shown in Figure [l2| (right; r is the correlation coefficient). Clusters with CC, WCC, and NCC 
are plotted with blue filled circles, green open circles, and red open triangles, respectively. There is a clear anti-correlation between the 
Lx~r s i scatter and the cooling time, so that CC clusters tend to lie above the mean relation. 



i Sijacki et al.||2007 l, and to inject mass and momentum as 
well in the form of bipolar jets ( |Cattaneo fc Teyssier|2007l ). 
For the survey of parameter sensitivity, we investigated pa- 
rameters that are least constrained by observations, includ- 
ing both numerically relevant parameters (Table[T]and Table 
[2| and physically motivated parameters (Table |3^. We com- 
pared their influence on the evolution of SMBH and cluster 
properties and examined their ability to self-regulate and 
reproduce observed profiles inside cluster cores. The main 
findings of the sensitivity study are summarized in the fol- 
lowing. 

1. Resolution - The convergence tests show that in- 
creasing resolution generally produces more variable accre- 
tion rates. The bubble feedback suffers greater variation by 
changing the resolution, whereas the jet model is more ro- 



bust, as long as the radius for computing accretion rates is 
larger than the sizes of the jets. 

2. Accretion - The proportionality used to relate the 
Bondi accretion to the actual SMBH accretion rate has a 
significant impact on the evolution of SMBH and cluster 
properties. Given the uncertainties in the accretion mecha- 
nisms, current AGN subgrid models may have very limited 
power to predict the evolution of cluster core properties, 
such as the fraction of CC versus NCC clusters as a func- 
tion of time. 

3. Efficiency of mechanical heating - Varying the me- 
chanical heating efficiency does not alter the overall evolu- 
tion as much as accretion. Feedback with large efficiencies 
has more variable accretion rates, more suppression in black 
hole growth, and higher entropy floors. Efficiencies that are 
too small would fail to overcome cooling. 
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4. Frequency - Changing the frequencies of injections 
has a minor effect. Longer duty cycles tend to generate more 
fluctuations in the accretion rates and cluster profiles. 

5. Region - The evolution of the SMBH and the ICM is 
very sensitive to the size of the energy injection region (the 
displacement from the BH does not matter much). More- 
over, bubbles that are too large would sometimes produce 
entropy floors that are inconsistent with observations, and 
bubbles that are too small may not be able to heat the entire 
CC and stop catastrophic cooling. Thus for any model that 
requires setting the feedback sizes by hand, there would be 
an issue of fine-tuning for a general population of clusters. 

6. Kinetic feedback - The jets with varied thermal to 
kinetic ratios produce very similar results. A combination 
of thermal and kinetic energy is slightly more efficient than 
purely thermal feedback. Purely kinetic feedback with effi- 
ciencies that are too small would fail to self-regulate. 

Comparing the bubble and jet models, we find that their 
main difference lies in the sizes of energy injection regions 
(Figure j8J|. The models are numerically degenerate when ap- 
propriate parameters are chosen, i.e., producing tiny, contin- 
uous bubbles to mimic the jets. The jet model is in general 
more robust to many numerical parameters (e.g. resolution) 
as well as physical parameters (e.g. sizes of feedback, which 
is desirable because there is no need for fine-tuning). How- 
ever, though the jets can maintain the cluster in the CC 
state, to avoid the artificial accumulation of cold gas around 
the black hole requires better treatment of the multiphase 
gas. Also, purely thermal concentrated heating, like thermal 
jets or quasar feedback, would produce a central peak in 
the temperature profile. Therefore, one needs to be cautious 
when interpreting results in the immediate surroundings of 
the BH. 

Outbursts from AGN are energetic events that can 
greatly influence the observable properties of galaxy clus- 
ters. Previous simulations with AGN subgrid models have ei- 



ther studied their impact inside cluster cores ( Gaspari et al 



2011 Dubois et al. 20101, or focused on matching the ob- 



served scalings of SMBH evolution and cluster gross prop- 



erties (jSijacki et al-pOCTf] |Booth fc Sch ayc 2009; |Puchwein 



et al. 12008 1. However, whether these models can simultane- 
ously reproduce cluster properties both inside and outside 
the cores has not previously been demonstrated. In the sen- 
sitivity study we identified model parameters that can self- 
regulate and produce core profiles consistent with observa- 
tions. Now we summarize our findings for cluster integrated 
properties as follows. 

1. All the subgrid models that successfully regulate cool- 
ing in the previous analysis also produce variation in inte- 
grated quantities consistent with the scatter of the observed 
scaling relations (Figure |9j. 

2. The model uncertainties in M s , T B \, Lx, Ysz, and Yx 
as functions of overdensity radius are quantified in Table [4] 
Quantities that are more sensitive to gas density (e.g. M s , 
Lx) have larger uncertainties, whereas T s i, Ysz, and Yx are 
most robust to model variations, to the levels of ~ 10 — 20% 
at 7?5oo, and ~ 5 — 10% at 7?200- 

3. Since AGN feedback reduces gas density and raises 
temperature, anti-correlations exist between Lx and T s i and 
also between M g and T B \, contributing to the intrinsic scatter 
in these two scaling relations. However, because the ICM 



reacts to AGN feedback with a delay, correlations between 
observables and AGN power are weak. 

4. Because M g and T s i are anti-correlated, even under 
the influence of strong AGN outbursts, the Ysz and Yx 
parameters are still robust mass proxies. 

Contrasting the bubble and jet models, we find that the 
more extended bubble injections are generally more effective 
in altering cluster properties than the more concentrated 
jet feedback. Consequently, simulations using the bubble 
feedback model and the accretion strength a = 100 ( |Puch-| 
wein et aL]|2008| are able to steepen the Lx~Tx slope but 
have difficulties producing CC clusters as observed. Though 
studies based on an improved accretion model (i.e. the ft 



model proposed by Booth & Schaye (20091) have success- 



fully matched the Lx~Tx slope and other properties on 
group scales (McCarthy et aL||2010[ ), future simulations on 
the cluster scale are required to verify whether CC clus- 
ters can be produced. If not, then it could mean that the 
sizes of the bubbles are still too large and have to be fur- 
ther controlled, or that other accretion models need to be 
considered. On the other hand, though simulations adopt- 
ing the jet model (either for an idealized cluster as in our 
study, or re-simulations from a cosmological volume as in 



Dubois et al. (20101 and Gaspari et al. (2011 1) can success- 
fully maintain clusters in the CC state, a statistical sample 
of clusters generated from a full cosmological simulation us- 
ing the jet model does not yet exist to verify whether the 
jets could provide enough entropy to steepen the Lx~Tx 
relation. If not, either the feedback energy needs to be dis- 
tributed by some other mechanism, or the solution still lies 
in other accretion models. We recommend that these possi- 
bilities should be investigated to understand the limitations 
of the existing bubble and jet models before one attempts 
to refine the parameter space of any particular model. 

The integrated Compton y parameter, Ysz, and its X- 
ray analog, Yx, are considered very good cluster mass prox- 
ies because previous simulations (without AGN feedback) 



show that they have very small mass scatter ( Kravtsov et al 



2006 1 and they are relatively insensitive to cluster dynamical 



state (Poole et aLpO"07| |Wik et al.|2008) |Yang et al-pOlo] ) 
Here we further show that the Y parameters are not easily 
disturbed by powerful AGN outbursts, which adds another 
reason to why they present so little observed scatter and can 
be used as excellent mass tracers. 

As it becomes more common to use the scaling rela- 
tions of the Y parameters provided by numerical simula- 



tions for calibration in observational studies (Arnaud et al 



20101 or for deriving cosmological constraints (Mantz et al 



20fJ8||Vikhlinin et al.|2009)|Vanderlinde et al.|2010| >, we note 
that incomplete knowledge of the processes of AGN feedback 
puts a limit on the predictive power of current cosmological 
simulations. Even for these most robust variables, the theo- 
retical uncertainties due to model variations are ~ 10 — 20% 
at -R500 and ~ 5 — 10% at -R200, which would translate into 
mass errors comparable to other main sources of systematic 
errors reported in the literature, such as the bias of hydro- 



static mass due to non-thermal pressure support (e.g. Lau 



|et al.|[2009[ ). Those variables that are sensitive to gas den- 
sity (e.g. Lx in particular) are even more uncertain. We note 
that the level of model uncertainties may be dependent on 
cluster masses, as the gas would be more easily displaced 
by AGN feedback in lower mass systems. For clusters more 
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massive than our simulated cluster, the model uncertain- 
ties of integrated quantities may be smaller than the values 
quoted above. However, since the predictions get progres- 
sively worse near cluster cores, it is very likely that there are 
still substantial uncertainties inside the cores of more mas- 
sive clusters. Furthermore, since clusters form hierarchically, 
the core properties (e.g. entropy floor) of massive clusters at 
the present day are determined by when and how the gas is 
heated inside the lower-mass systems that later merge into 
the clusters (e.g. McCarthy et al.|2011 1. Therefore, in order 
to improve predictions of cluster observables (including the 
cores) and derive cosmological constraints to the percent 
level, it is essential for numerical simulations to focus on 
how to improve the modeling of the subgrid physics before 
making various predictions. 

Our sensitivity study showed that in order to effec- 
tively improve future AGN subgrid models, the crucial next 
step is to further constrain the accretion processes of the 
SMBH, the mechanical heating efficiency, and the sizes of the 
feedback region. The mechanical heating efficiency param- 
eter (i.e. the fraction of feedback energy transformed into 
heat) and the sizes of the region to distribute heat may be 



evaluated from detailed numerical simulations (Vernaleo & 



Reynolds 


2007||O'Neill & Jones|2010|| 


in § |3.3.2 


the results would depend 



the simulations. Thus to pin down these parameters still re- 
quires more knowledge of the mixing properties of the ICM. 
For the mechanical heating efficiency, before its value can 
be estimated reliably, an alternative way is to adjust the ef- 
ficiency to match the normalization of the Mbh-c relation 
or the cosmic black hole mass density at z = ( Sijacki et al. 
2007||Booth fe Schaye|2009| |Dubois et al.|2011| >. Note that 
since the obtained value is also dependent upon the specific 
accretion and feedback model employed, the efficiency pa- 
rameter, if determined in this fashion, cannot inherit from 
other simulations but will have to be normalized for each 
realization. 

Improving the subgrid accretion model is a more chal- 
lenging task, simply because it is still unclear how to link the 
accretion rates across such a great dynamical range. So far 
most cosmological simulations evaluate accretion rates based 
on the Bondi accretion rate. While the Bondi accretion rate 
appears sufficient in powering the observed AGN jets for 
many cases ( | Allen et al.|2006| ), for some systems other mech- 



anisms seem to be required (McNamara et al. 20111. Fur- 



thermore, the original Bondi accretion rate is based on sim- 
plified assumptions such as spherically-symmetric, steady 
flow with zero velocity at the Bondi radius. These crite- 
ria may not be applicable to all systems, especially those 
with radial infall due to rapid cooling, or with non-negligible 
angular momentum (e.g. Power et al. 20111. Alternative 



schemes have been proposed, including stochastic accretion 
( |Pope||2007| ), cold gas accretion ( jPizzolato fc Soker||2005[), 
accretion by gravitational instabilities in galaxies (Hopkins 
|& Quataert|2010l), and SMBH spins dMcNamara et al. 120091). 



These models are not yet integrated into cosmological simu- 



lations (except a cold-accretion-like scheme used in Gaspari 



et al. (2011 1). Clearly more detailed investigations and com- 
parisons in this area are necessary for further improvement 
of the AGN subgrid models. 
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